R Notebook

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages ------------------------------------------------------------------- tidyverse 1.3.1 --
v ggplot2 3.3.5     v purrr   0.3.4
v tibble  3.1.2     v dplyr   1.0.7
v tidyr   1.1.3     v stringr 1.4.0
v readr   1.4.0     v forcats 0.5.1
-- Conflicts ---------------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(spatstat)

Lade nötiges Paket: spatstat.data
Lade nötiges Paket: spatstat.geom
Registered S3 method overwritten by 'spatstat.geom':
  method     from
  print.boxx cli 
spatstat.geom 2.3-0
Lade nötiges Paket: spatstat.core
Lade nötiges Paket: nlme

Attache Paket: ‘nlme’

Das folgende Objekt ist maskiert ‘package:dplyr’:

    collapse

Lade nötiges Paket: rpart
spatstat.core 2.3-0
Lade nötiges Paket: spatstat.linnet
spatstat.linnet 2.3-0

spatstat 2.2-0       (nickname: ‘That's not important right now’) 
For an introduction to spatstat, type ‘beginner’ 

library(pracma)

Attache Paket: ‘pracma’

Das folgende Objekt ist maskiert ‘package:spatstat.core’:

    rat

Die folgenden Objekte sind maskiert von ‘package:spatstat.geom’:

    integral, quad

Das folgende Objekt ist maskiert ‘package:purrr’:

    cross

library(plyr)

----------------------------------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------------------------------------------------------------

Attache Paket: ‘plyr’

Die folgenden Objekte sind maskiert von ‘package:dplyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

Das folgende Objekt ist maskiert ‘package:purrr’:

    compact

library(plotly)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attache Paket: ‘plotly’

Die folgenden Objekte sind maskiert von ‘package:plyr’:

    arrange, mutate, rename, summarise

Das folgende Objekt ist maskiert ‘package:ggplot2’:

    last_plot

Das folgende Objekt ist maskiert ‘package:stats’:

    filter

Das folgende Objekt ist maskiert ‘package:graphics’:

    layout

library(mltools)

Attache Paket: ‘mltools’

Das folgende Objekt ist maskiert ‘package:tidyr’:

    replace_na

library(data.table)

data.table 1.14.0 using 2 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attache Paket: ‘data.table’

Das folgende Objekt ist maskiert ‘package:spatstat.geom’:

    shift

Die folgenden Objekte sind maskiert von ‘package:dplyr’:

    between, first, last

Das folgende Objekt ist maskiert ‘package:purrr’:

    transpose

library(htmlwidgets)

Function to calculate B1 and B2 of a specific omega.

calcB2 <- function(X1w, X2w, X3w) {
  norm12 = norm(X1w-X2w,"2")
  norm13 = norm(X1w-X3w, "2")
  norm23 = norm(X2w-X3w, "2")
  N = c(norm12, norm13, norm23)
  ind <- which.max(N)
  if(ind == 1){
    B2w = X3w
  }else if(ind == 2){
    B2w = X2w
  }
  else{
    B2w = X1w
  }
  return(list("B2w" = B2w, "ind" = ind))
}

Calculate ordered random variables for a specific omega.

calcOrder <- function(X1w, X2w, X3w, B2w, ind){
  norm1 = norm(X1w-B2w, "2")
  norm2 = norm(X2w-B2w, "2")
  norm3 = norm(X3w-B2w, "2")
  X3wP <- B2w
  if (ind == 1){
    if(norm2 >= norm3){
      X1wP <- X2w
      X2wP <- X3w
    }else{
      X1wP <- X3w
      X2wP <- X2w
    }
  }else if(ind == 2){
    if(norm1 >= norm3){
      X1wP <- X1w
      X2wP <- X3w
    }else{
      X1wP <- X3w
      X2wP <- X1w
    }
  }else{
    if(norm1 >= norm2){
      X1wP <- X1w
      X2wP <- X2w
    }else{
      X1wP <- X2w
      X2wP <- X1w
    }
  }
  
  return(list("X1" = X1wP, "X2" = X2wP, "X3" = X3wP))
}

Calculate F_i for a specific Omega

calcF <- function(R2, B2w, X1wP, X2wP, X3wP){
  
  # calculate F1
  # browser()
  n <- length(R2$x)
  F1ind = numeric()
  b <- Norm(X1wP - B2w)
  for(i in 1:n){
    # browser()
    x = as.numeric(R2[i,])
    a <- Norm(x - B2w)
    if(a >= b){
      F1ind <- c(F1ind,i)
    }
  }
  F1 <- R2[F1ind,]
  #R2 <- R2[-F1ind,]
  
  # calculate F2
  F2ind = numeric()
  b <- Norm(X2wP - B2w)
  # print(b)
  if(n!=0){
    for(i in 1:n){
      if(!(i %in% F1ind)){
        x = as.numeric(R2[i,])
        a <- Norm(x - B2w)
        if(a >= b){
          F2ind <- c(F2ind,i)
        }
      }
    }
  }
  F2 <- R2[F2ind,]
  #R2 <- R2[-F2ind,]
  
  #calculate F3
  #n <- length(R2$x)
  F3ind = numeric()
  b <- Norm(X3wP - B2w)
  # print(b)
  if(n != 0){
    for(i in 1:n){
      if(!(i %in% F1ind) && !(i %in% F2ind)){
        x = as.numeric(R2[i,])
        a <- Norm(x - B2w)
        if(a >= b){
          F3ind <- c(F3ind,i)
        }
      }
    }
  }
  F3 <- R2[F3ind,]
  #R2 <- R2[-F3ind,]
  
  #calculate F4 (the rest) and F0 (the empty set) 
  F4 <- R2
  
  return(list("F1ind" = F1ind, "F2ind" = F2ind, "F3ind" = F3ind))
}

Calculate Probability of one i which is part of C

calcProb <- function(sizeOmega, ind_x, ind, allF){
  count = 0
  for(i in 1:sizeOmega){
    Fi_ind <- allF[ind,i]
    # print(x)
    # print(nrow(merge(x,Fi)))
    Fi_ind = unlist(Fi_ind)
    if(ind_x %in% Fi_ind){
      count = count + 1
    }
  }
  # print(count/sizeOmega)
  return(count/sizeOmega)
}

Calculate C for one x (F_i (x) still missing)

funcC <- function(ind_x, sizeOmega, allF, empCDF){
  res = 0
  
  for(ind in 1:3){
    # print(res)
    res = res + calcProb(sizeOmega, ind_x, ind, allF)*(empCDF[ind_x]+(ind-1))
  }
  # print(res)
  return(res/4)
}

calcAllC <- function(){
  allF = replicate(sizeOmega, calcAllF())
  n <- length(R2$x)
  Cxy <- numeric(n)
  dt <- data.table(x = rnorm(100,-5,5), y = rnorm(100,-5,5))
  empCDF <- empirical_cdf(dt, ubounds=CJ(x = xR2, y = yR2))
  empCDF <- empCDF$CDF
  for(i in 1:n){
    Cxy[i] <- funcC(i, sizeOmega, allF, empCDF)
  }
  final <- R2
  final$Cxy <- Cxy
  return(final)
}

Calculate all the X_(i), F_i for all omega

calcAllF <- function(){
  # X1w = runif(2,lowerBound, upperBound)
  # X2w = runif(2,lowerBound, upperBound)
  # X3w = runif(2,lowerBound, upperBound)
  X1w = rpois(2,4)
  X2w = rpois(2,4)
  X3w = rpois(2,4)
  
  ret <- calcB2(X1w,X2w,X3w)
  B2w <- ret$B2w
  ind <- ret$ind
  
  ret <- calcOrder(X1w, X2w, X3w, B2w, ind)
  X1wP <- ret$X1
  X2wP <- ret$X2
  X3wP <- ret$X3
  
  scriptF <- calcF(R2, B2w, X1wP, X2wP, X3wP)
  
  return (scriptF)
}

Some testing